In phuonghtht/PMA_dengue_v2: Lab diagnosis

loading saved fits

# Ref:https://cran.r-project.org/web/packages/loo/vignettes/loo2-with-rstan.html
#https://www.jmp.com/support/help/en/15.2/index.shtml#page/jmp/likelihood-aicc-and-bic.shtml
#https://mc-stan.org/rstan/reference/stan.html

rm(list=ls())
library(loo)
library(rstan)

# loading the model fit_MLE
site = "HC";dataInput = "C1D1"
model <- "TVSS" # TV, SS, TVSS
groupAge <-1 # HC:137  KH:125
step_year <- 1
dataYear <- 35
if(site =="HC") {
  dataPoint <- 137
}
if(site =="KH") {
  dataPoint <- 125
}

if(model =="cons"|model == "SS"){
  fit <- readRDS(paste("/home/phuong/pCloudDrive/PMA analysis results/FoI/RDS/FOI",model,"groupAge",groupAge,site,dataInput, '.rds'))
}
if(model =="TV"|model == "TVSS"){
  fit <- readRDS(paste("/home/phuong/pCloudDrive/PMA analysis results/FoI/RDS/FOI",model,"groupAge",groupAge,"step_year",step_year,site,dataInput, ".rds"))
}


# fit1<- readRDS(paste("/home/phuong/pCloudDrive/PMA analysis results/FoI/RDS/FOI cons groupAge",groupAge,site, dataInput, '.rds'))
# fit2 <- readRDS(paste("/home/phuong/pCloudDrive/PMA analysis results/FoI/RDS/FOI  TV groupAge",groupAge,"step_year",step_year,site, dataInput, '.rds'))
# fit3<- readRDS(paste("/home/phuong/pCloudDrive/PMA analysis results/FoI/RDS/FOI  SS groupAge",groupAge,site, dataInput, '.rds'))
# fit4<- readRDS(paste("/home/phuong/pCloudDrive/PMA analysis results/FoI/RDS/FOI TVSS groupAge",groupAge,"step_year",step_year,site, dataInput, '.rds'))


# Extract pointwise log-likelihood
# using merge_chains=FALSE returns an array, which is easier to 
# use with relative_eff()

loglike = extract_log_lik(fit, parameter_name = "log_like", merge_chains = T)

# loglike1 = extract_log_lik(fit1, parameter_name = "log_like", merge_chains = T)
# loglike2 = extract_log_lik(fit2, parameter_name = "log_like", merge_chains = T)
# loglike3 = extract_log_lik(fit3, parameter_name = "log_like", merge_chains = T)
# loglike4 = extract_log_lik(fit4, parameter_name = "log_like", merge_chains = T)

Extracting data from the fits

# fit <- fit1
# util$check_all_diagnostics(fit)

posterior <- rstan::extract(fit) # specify package rstan rather than tidyr

# posterior1 <- rstan::extract(fit1) # specify package rstan rather than tidyr
# posterior2 <- rstan::extract(fit2) # specify package rstan rather than tidyr
# posterior3 <- rstan::extract(fit3) # specify package rstan rather than tidyr
# posterior4 <- rstan::extract(fit4) # specify package rstan rather than tidyr

# plot(posterior1$lambda,posterior1$sumloglike)


medianLL <- median(posterior$sumloglike)

# medianLL1 <- median(posterior1$sumloglike)
# medianLL2 <- median(posterior2$sumloglike)
# medianLL3 <- median(posterior3$sumloglike)
# medianLL4 <- median(posterior4$sumloglike)

# meanLL1 <- mean(posterior1$sumloglike)
# meanLL2 <- mean(posterior2$sumloglike)
# meanLL3 <- mean(posterior3$sumloglike)
# meanLL4 <- mean(posterior4$sumloglike)

Comparing models

# function for calculating AIC
cal.AIC.c <- function(LL,k,n){
  #LL: loglikelihood
  #k:the number of estimated parameters in the model
  #n: the number of observations used in the model
  AIC = -2*LL +2*k + 2*k*(k+1)/(n - k -1)
  print(AIC)
}
cal.AIC <- function(LL, k){
  AIC = -2*LL + 2*k
  print(AIC)
}

if(model == "cons"){
  AIC <- cal.AIC( LL= median(posterior$sumloglike), k = 1)
  AICc <- cal.AIC.c( LL= median(posterior$sumloglike), k = 1, n=dataPoint)
}
if(model == "TV"){
  AIC2 <- cal.AIC( LL= median(posterior$sumloglike), k = ceiling(dataYear/step_year))
  AICc2 <- cal.AIC.c( LL= median(posterior$sumloglike), k = dataYear/step_year, n=dataPoint)#TV 
}
if(model == "SS"){
  AIC <- cal.AIC( LL= median(posterior$sumloglike), k = 4)
  AICc <- cal.AIC.c( LL= median(posterior$sumloglike), k = 4, n=dataPoint)
}
if(model == "TVSS"){
  AIC <- cal.AIC( LL= median(posterior$sumloglike), k = ceiling(dataYear/step_year)*4)#TVSS
  AICc <- cal.AIC.c( LL= median(posterior$sumloglike), k = ceiling(dataYear/step_year)*4, n=dataPoint)#TVSS
}


# ##AIC :
# AIC1 <- cal.AIC( LL= median(posterior1$sumloglike), k = 1)#cons
# AIC2 <- cal.AIC( LL= median(posterior2$sumloglike), k = ceiling(dataYear/step_year))#TV
# AIC3 <- cal.AIC( LL= median(posterior3$sumloglike), k = 4)#SS
# AIC4 <- cal.AIC( LL= median(posterior4$sumloglike), k = ceiling(dataYear/step_year)*4)#TVSS
# ##AICc :
# AICc1 <- cal.AIC.c( LL= median(posterior1$sumloglike), k = 1, n=dataPoint)#cons
# 
# AICc2 <- cal.AIC.c( LL= median(posterior2$sumloglike), k = dataYear/step_year, n=dataPoint)#TV
# AICc3 <- cal.AIC.c( LL= median(posterior3$sumloglike), k = 4, n=dataPoint)#SS
# 
# AICc4 <- cal.AIC.c( LL= median(posterior4$sumloglike), k = ceiling(dataYear/step_year)*4, n=dataPoint)#TVSS